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Abstract 

In this paper, we study decoherence in Grover's quantum search algorithm using 
a perturbative method. We assume that each two-state system (qubit) that belongs 
to a register suffers a phase flip error (a z error) with probability p independently 
at every step in the algorithm, where < p < 1. Considering an n-qubit density 
operator to which Grover's iterative operation is applied M times, we expand it 
in powers of IMnp and derive its matrix element order by order under the large-n 
limit. [In this large-n limit, we assume p is small enough, so that 2Mnp can take any 
real positive value or zero. We regard x = 2Mnp(> 0) as a perturbative parameter.] 
We obtain recurrence relations between terms in the perturbative expansion. By 
these relations, we compute higher orders of the perturbation efficiently, so that 
we extend the range of the perturbative parameter that provides a reliable analysis. 
Calculating the matrix element numerically by this method, we derive the maximum 
value of the perturbative parameter x at which the algorithm finds a correct item 
with a given threshold of probability P t h or more. (We refer to this maximum value 
of x as x c , a critical point of x.) We obtain a curve of x c as a function of Pth by 
repeating this numerical calculation for many points of Pth and find the following 
facts: a tangent of the obtained curve at P t h = 1 is given by x = (8/5)(l — Pth), 
and we have x c > —(8/5) log e P t h near P t h = 0. 



1 Introduction 

Many researchers think that decoherence is one of the most serious difficulties in realizing 
the quantum computation [H El EH E] • The decoherence is caused by interaction between 
the quantum computer and its environment. The interaction lets the state of the computer 
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become correlated with the state of the environment. Consequently, some of information of 
the quantum computer leaks into the environment. This process causes errors in the state 
of the quantum computer, and as a result, the probability that the quantum algorithm 
gives the right answer decreases. To overcome this problem, quantum error-correcting 
codes are proposed 0E1IZI- 

Not only for practical purposes but also for theoretical interests, an important question 
is how robust the quantum algorithm is against this disturbance. If we know the upper 
bound of the error rate that allows the quantum computer to obtain a solution with a 
certain probability or more, this bound is useful for us to design quantum gates. 

Grover's algorithm is considered to be an efficient amplitude-amplification process for 
quantum states. Thus it is often called a search algorithm jHlEl- By applying the same 
unitary transformation to the state in iteration and gradually amplifying an amplitude of 
one basis vector that an oracle indicates, Grover's algorithm picks it up from a uniform 
superposition of 2 n basis vectors with a certain probability in 0(2 n//2 ) steps. In view 
of computational time (the number of queries for the oracle), the efficiency of Grover's 
algorithm is proved to be optimal |lUj . 

In Ref. jTTj, we study decoherence in Grover's algorithm with a perturbative method. 
We consider the following simple model. First, we assume that we search |0...0) from 
the uniform superposition of 2 n logical basis vectors {\x) : x G {0,l} n } by Grover's 
algorithm. This assumption simplifies the iterative transformation. Second, we assume 
that each qubit of the register interacts with the environment independently and suffers 
a phase damping, which causes a phase flip error (a z error) with probability p and does 
nothing with probability (1— p) to the qubit. In this model, we expand an n-qubit density 
operator to which Grover's iterative operation is applied M times in powers of 2Mnp. 
Then, we take the large-n limit, so that we can simplify each order term of the expansion 
of the density operator and we obtain its asymptotic form. 

In this large-n limit, we assume p is small enough, so that 2Mnp can take any real 
positive value or zero. We regard x = 2Mnp(> 0) as a perturbative parameter. We can 
interpret x = 2Mnp as the expected number of phase flip errors (a z errors) that occur 
during the running time of computation. In Ref. [TT], we give a formula for deriving an 
asymptotic form of an arbitrary order term of the perturbative expansion. However, this 
formula includes a complicated multiple integral and the number of terms in its integrand 
increases exponentially. Because of these difficulties, we obtain explicit asymptotic forms 
only up to the fifth-order term. 

In this paper, using recurrence relations between terms of the perturbative expansion, 
we develop a method for computing higher order terms efficiently. By this method, we 
derive an explicit form of the density matrix of the disturbed quantum computer up to 
the 39th-order term with the help of a computer algebra system. (In actual fact, we 
use Mathematica for this derivation.) Because we consider the higher order perturbation, 
we can greatly extend the range of the perturbative parameter that provides a reliable 
analysis, compared with our previous work in Ref. [TT]. Calculating the matrix element 
up to the 39th-order term numerically from the form obtained by the computer algebra 
system, we derive the maximum value of the perturbative parameter x at which the 
algorithm finds a correct item with a given threshold of probability P t ^ or more. (We 
refer to this maximum value of x as i c , a critical point of x.) 
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Figure 1: A schematic representation of function of P t h- -Pth is a threshold of 

probability, x represents both the perturbative parameter and the expected number of 
errors during the running time of computation. x c is a critical point of x. Both P t h and 
x are dimensionless. We can easily obtain x c = for P t h = 1- This fact is included in the 
above schematic graph. The above graph represents a phase diagram that consists of two 
domains. One domain is where the quantum algorithm is effective and the other domain 
is where it is not effective. 

Grover's algorithm can find the correct item by less than (7r/4)v2™ steps with given 
probability P^ or more under no decoherence (p = 0). When we fix P t h, the number of 
iterations that we need increases as the decoherence becomes stronger (p becomes larger). 
Finally we never detect the correct item with P t h or more for p > p c . Thus, we can 
think p c to be a critical point for P th . (p c depends on P t h-) However, we actually obtain 
x c = x c(Pth) for the perturbative parameter x = 2Mnp instead of p c = p c (P t h)- From 
the relation x c = x c (P t h), we can draw a phase diagram as shown in Fig. ^ The diagram 
consists of two domains. One is where the quantum algorithm is effective and the other 
is where it is not effective. 

Figures El and El represent a curve of x c = x c (Pth) obtained by repeating the numerical 
calculation of x c for many points of P t h- In Fig. (3, we use a linear scale on both horizontal 
and vertical axes. We prove later that a tangent of the curve x = x c (P t h) at P t h = 1 is 
given by x = (8/5)(l — Pth)- In Fig. El we use log and linear scales on the horizontal and 
vertical axes, respectively. We observe x c > —(8/5) log e P t h nea r P t h = from this figure. 

Here, we mention that we can investigate our model by Monte Carlo simulations, 
as well. In fact, we compare results obtained by our perturbative method with results 
obtained by Monte Carlo simulations in Figs. Eland HI in Sec. 01 and we confirm that they 
are consistent. From these analyses, we conclude that our perturbative method is valid 
in a certain range of the perturbative parameter. 

However, the Monte Carlo simulation method has some difficulties for investigating our 
model. First, the execution time of computation increases exponentially in n (the number 
of qubits). We have to always come up against this problem when we simulate a process 
of a quantum computer with a classical computer. Second, the Monte Carlo simulation 
method is not suitable for obtaining a variation of a physical quantity as a function of 
some parameters, because we carry out each simulation with fixing parameters such as 
the error rate p and the threshold of probability P t h- Thus we prefer our perturbative 
method to the Monte Carlo simulation method for computing x c (the critical point of x) 
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Figure 2: x c as a function of P t h- [A thick solid curve represents x = x c (P t h)-] -Pth is a 
threshold of probability. x c is a critical point of x (the perturbative parameter). We use 
a linear scale on both horizontal and vertical axes. The data are obtained by repeating 
numerical calculation of x c for many points of P t h- Because x = x c (P t h) shows a sharp 
divergence at P t h = 0, we start calculation of x c from P t h = 1. While we are going from 
P t h = 1 toward P t h = 0, we make a finite difference of P t h smaller gradually. (We put 
AP th = 5.0 x 10" 4 around P th = 1 and AP th = 5.0 x 10~ 7 around P th = 3.7 x 10~ 3 .) A 
thin dashed line represents a tangent of x = x c (P t h) at P t h = 1- 




Figure 3: x c as a function of P t h- [A thick solid curve represents x = x c (P t h).] We use log 
and linear scales on the horizontal and vertical axes, respectively. In this figure, we use the 
same data of x = x c (P t h) as in Fig. |21 A thin dashed line represents x = —(8/5) log e P t h- 
We observe x c > —(8/5) log e P t h near P t h = 0. 
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that is obtained by evaluating the probability of detecting a correct answer as a function 
of x and P t h- 

A related result is obtained in the study of the accuracy of quantum gates by Bernstein 
and Vazirani jT^j, and Preskill They consider a quantum circuit where each quantum 
gate has a constant error because of inaccuracy. Thus, it is an error of a unitary trans- 
formation and it never causes dissipation of information from the quantum computer to 
its environment. They estimate inaccuracy e for which the quantum algorithm is effective 
under the fixed number of time steps T, and obtain 2Te < 1 — P t h> where < e < 1. If 
we regard p/2 as inaccuracy e and 2Mn as the number of whole steps in the algorithm 
T, it is similar to our observation that x c = 2Mnp ~ (8/5) (1 — P t h) near P t h = 1, except 
for a factor. 

Barenco et al. study the approximate quantum Fourier transformation (AQFT) and 
its decoherence jHj. Although their motivation is slightly different from Refs. ^21 Eh we 
can think their model to be the quantum Fourier transformation (QFT) with inaccurate 
gates. They confirm that AQFT can give a performance that is not much worse than the 
QFT. 

This article is organized as follows. In Sec. EJ we describe our model and perturbation 
theory defined in our previous work jTTj. In Sec. El we give recurrence relations between 
terms of the perturbative expansion. We develop a method for calculating higher order 
perturbation efficiently with these relations. In Sec. EJ we carry out numerical calculations 
of the matrix element of the density operator by the efficient method obtained in Sec. El 
Moreover, we investigate the critical point x c and obtain the phase diagram shown in 
Figs. 121 and El In Sec. El we give brief discussions. In Appendix^ we give a proof of an 
equation which appears in Sec. El 



2 Model and perturbation theory 

In this section, we first describe a model that we analyze. It is a quantum process of 
Grover's algorithm under a phase damping at every iteration. Second, we formulate a 
perturbation theory for this model. 



2.1 Model 

First of all, we give a brief review of Grover's algorithm ,8j. Starting from the ra-qubit 
uniform superposition of logical basis vectors, 

W\Q...0) = -^= \ x ) forn>2, (1) 

v 2 n £(={0,1}™ 

Grover's algorithm gradually amplifies the amplitude of a certain basis vector \xq) that 
a quantum oracle indicates, where xq G {0, 1}". The operator W in Eq. is an n- 
fold tensor product of a one-qubit unitary transformation and given by W = H® n . The 
operator H is called Hadamard transformation and represented by the following matrix, 



» = ^ _\ , (2) 




5 



where we use the orthonormal basis {|0), |1)} for this matrix representation. The quantum 
oracle can be regarded as a black box and actually it is a quantum gate that shifts phases 
of logical basis vectors as 

R ■ I ^ ~ * ~\ Xo ^ (3) 

X ° \ \ x ) ~^ \ x ) f° r x 7^ ^0 ' 

where xq,x G {0, l} n . [We note that all operators (quantum gates) in Grover's algorithm 
are unitary. Thus, W = H" 1 , = W~ l , R\, Q = R^, and so on.] 

To let probability of observing \xq) be greater than a certain value (1/2, for example), 
we repeat the following procedure 0(v2") times: 

1. Apply R Xo to the n-qubit state. 

2. Apply D = WR W to the n-qubit state. 

Ro is a selective phase-shift operator, which multiplies a factor (—1) to |0...0) and does 
nothing to the other basis vectors, as defined in Eq. (J5|i. D is called the inversion-about- 
average operation. 

From now on, we assume that we amplify an amplitude of |0...0). From this assump- 
tion, we can write an operation iterated in the algorithm as 

DR = (WR W)R . (4) 

After repeating this operation M times from the initial state W|0)(= VV|0...0)), we obtain 
the state (WR ) 2M W\0). (We often write |0) as an abbreviation of the n-qubit state |0...0) 
for a simple notation.) 

Next, we think about the decoherence. In this paper, we consider the following one- 
qubit phase damping [13 EE] : 

P^ p' = V^zpOz + (1 - p)p for < p < 1, (5) 

where p is an arbitrary one-qubit density operator. a z is one of the Pauli matrices and 
given by 

where we use the orthonormal basis (|0), |1)} for this matrix representation. For sim- 
plicity, we assume that the phase damping of Eq. © occurs in each qubit of the register 
independently before every Rq operation during the algorithm. This implies that each 
qubit interacts with its own environment independently. 

Here, we add some notes. First, because Ro G U{2 n ) is applied to all n qubits and 
H G U(2) is applied to only one qubit, we can imagine that the realization of Rq is more 
difficult than that of W = H® n . Hence, we assume that the phase damping occurs only 
before Rq. Second, although we assume a very simple decoherence defined in Eq. (JSJ), we 
can think of other complicated disturbances. For example, we can consider decoherence 
caused by an interaction between the environment and two qubits and it may occur with 
a probability of 0(p 2 ). In this paper, we do not assume such complicated disturbances. 
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2.2 Perturbation theory 

Let p^ M > be the density operator obtained by applying Grover's iteration M times to the 
ate W 

We can expand p^ in powers of p and (1 — p) as follows: 



n-qubit initial state W|0). The decoherence defined in Eq. (JHJ) occurs 2Mn times in p^ M \ 



pM = (l-pf Mn T^ M) + {l-p) 2Mn - l P T[ M) + ... 

2Mn 

= E(l"P) 2Mn "VTi M) , (7) 

fc=0 



where } are given by 

T (M) = (WRo) 2M W\0){0\W(RoW) 2M , (8) 

n 2M-1 

T (M) = ^2 (WR ) 2M - l af\WRo) l W\0)(0\W(R W) l a { ;\R W) 2M - 1 , (9) 

i=l 1=0 
n n 2M-1 

T ( M) = J2J2J2(WR ) 2M - l a^\WR ) l W\0)(R.c.\ 

i=l 3=1 2=0 

i<3 

n n 2M-X2M-1-1 

+EE E E (^o) 2A/ -^ m cr«(^o) m ^ ) (^o)^|o) 

1=1 j = l 1=0 771=1 

x(H.c.|, (10) 

and so on, a-W represents the operator applied to the ith qubit for 1 < i < n, and (H.c. | 
represents a Hermitian conjugation of the ket vector on its left side. (Here, we note 
= W , Rq = R , and crW = a&).) We can regard as a density operator whose 
trace is not normalized. It represents the sum of states where k errors occur during the 
iteration of M operations. 

On the other hand, from Eq. (JJJ), we can expand p^ in powers of p as follows: 

p (M) = p (M) + 2Mnpp (M) + \ 2 Mn)(2Mn-l)p 2 pf I) + ... 

= E"f 2 f)rf", 2 (ii) 



s I * r 



where 



(M) _ rp(M) 
PO — J 5 

(M) _ T (M) £l 



k 



(2Mn\ 1 (k^ 



pl M) = (-l) fc E(-l)T forA; = 0,l,...,2Mn. (12) 



\ i / 



Here, let us take the limit of an infinite number of qubits (the large-n limit). We assume 
that we can take very small p, so that 2Mnp can be an arbitrary real positive value or 
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zero. If 2Mnp is small enough, we can consider x = 2Mnp(> 0) to be a perturbative 
parameter and the series of Eq. (fTTj) to be a perturbative expansion. 

Under this limit, we derive an asymptotic form of (0|p( A/ )|0). In the actual derivation, 
we take the limit of n — ► oo with holding x = 2Mnp finite. (0|p^|0) is a probability 
that the quantum computer finds a correct item after M operations. Because we divide 
Tj by (Mn,y as in Eq. (JT2J), an expectation value of p^. can converge to a finite value 
in the limit n — > oo for k = 0,1, 2Mn. 

With these preparations, we will investigate the following physical quantities. Let P t h 
be a threshold of probability for < P t h < 1, so that if the quantum computer finds a 
correct item (in our model, it is |0)) with probability P t h or more, we regard it effective, 
and otherwise we do not consider it effective. Then, we consider the least number of the 
operations that we need to repeat for amplifying the probability of observing |0) to P t h 
or more for a given p. We refer to it as M th (p, P t h)- [M th (p, P th ) is the least number of M 
that satisfies (0|p (A/) |0) = P th for a given p.] As p becomes larger with fixing P t h, we can 
expect that M t h(p, Pth) increases monotonically. In the end, we never observe |0) at least 
with a probability P t h for a certain p c or more. (Hence, p c depends on Pth-) Regarding 
Pth as a threshold for whether the quantum computer is effective or not, we can consider 
p c to be a critical point. 

In our perturbation theory, we calculate physical quantities using the dimensionless 
perturbative parameter x = 2Mnp. Thus, we take M and x for independent variables. 
(In our original model defined in Sec. 12.11 we take M and p for independent variables.) We 
can define as well M t h(x, P t h) that represents the least number of the operations iterated 
for amplifying the probability of |0) to P t h for given x. Furthermore, we also obtain x c or 
more for which we can never detect |0) at least with probability P t h- 

Next, we evaluate (0|p ( - M ^|0). First, from simple calculation, we obtain the unper- 
turbed matrix element, 

(0|P (M) |0) P=0 = (0|T (M) |0) 

= sin 2 [(2M + 1)6], (13) 

where 



1 2 n - 1 

sin9 = ^, cos0 = W . 14) 

\/2P V 2 n 



(This parameter 9 is introduced by Boyer et al. [9 .) From Eq. 1)13)1. we notice the following 
facts. If there is no decoherence (p = 0), we can amplify the probability of observing |0) 
to unity. Taking large (but finite) n, we obtain sin9 ~ 9 and 9 ~ l/v2^, and we can 
observe |0) with unit probability after repeating Grover's operation M max ~ (7r/4)y / 2™ 
times. 

To describe the asymptotic forms of matrix elements, we introduce the following nota- 
tion. Because (0|T() |0) is a periodic function of M and its period is about tt\/2™ under 
the large-n limit, it is convenient for us to define a new variable G = lim^oo M9 (radian). 
Here, we give a formula for the asymptotic form of the kth order of the matrix element 
under n — > oo for k = 1,2, .... (The derivation of this formula is given in Sees. 4-7 and 
Appendix A of Ref. [Hj.) Preparing an fc-digit binary string a = (at, ...,ctk) £ {0, l} fe , 
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I T ai ,...,a k (<f>U ■■■Ak)\ 2 

r f sin 1 . , , I cos I I cos I ,„ , . I cos 



we define the following 2 k terms 

- n f (20i)rr a ) (2fa)..r:°\ w> u )\t\ we-T,*.))? 

cos sin sin sin ^ k ~i 

for k = 1,2,..., (15) 

where 



,e 








/ d 


'01 / 




/ 


/o 


JO 




JO 



and © denotes the addition modulo 2. We notice that the function of <j)\ and the other 
functions of fa,---, 4>k, @ — S«=i 0s are different (sine and cosine functions are put in 
reverse). These terms are integrated as 

®— 01-. --4>k-l 

d<fik 

x E l^,...,^,...,^)! 2 - (17) 

(Qi,...,a fc )e{0,l} fc 

We can obtain the matrix elements as follows. From Eq. (|13jh we obtain 

lim(0|T (M) |0) = sin 2 20. (18) 

n— >oo 

From Eqs. (fTKJ) and (fTTjl . we obtain 

/n|T( M )|n\ 1 r e 

lirn X ' ; = -j d0{[sin20cos2(6-0)] 2 + [cos 20 sin 2(6 - 0)] 2 } 

= 1-1^46-^^46, (19) 



(0|T 2 (M) |0) 
(Mn) 2 



lim 

n^oo 



1 

62 



1 

4 



e 



e-<t> 



dip 



o Jo 
x { [sin 20 cos 2y> cos 2(6 
+ [cos 20 cos2<£> sin2(6 - 
+ [sin20sin299sin2(6 — 
+ [cos 20 sin 2<p cos 2(6 - 



1 

16 



cos 46 



646 



sin 46, 



•-v?)] 2 
^)] 2 } 



(20) 



and so on. 
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The asymptotic form of the perturbative expansion of the whole density matrix is 
given by 

(P)(e,x) = lim(0|p (A/) |0) 

= C (Q) + C 1 (G)x+^C 2 (Q)x 2 + ... 

OO 1 

= E^)/ ( 21 ) 

fc=0 

where 

C (8) = F (6), 

d(e) = -F {e) + ^F 1 (e), 
k i A-i 

C*(e) = (-l) fc E(-2) J 7T— for fc = 0, 1, (22) 

F fc (e) = lim V ' * ' ; for fc = 0, 1, .... (23) 

In Eq. (J21|) . the fcth-order term is divided by fc!, so that we can expect the series (P)(0, x) 
to converge to a finite value for large x. 



and 



3 Recurrence relations between order terms 

In this section, we obtain recurrence relations between order terms of the perturbative 
series. Using this result, we develop a method for computing higher order terms efficiently 
When we compute lim^^^lT^^ |0) /(Mn) k for large k from Eqs. (fT5j) and (fTTj) . we 
notice the following difficulties: 

1. Equation (j!7|) includes an /cth-order integral. 

2. Equation (J 17)) includes 2 k terms being integrated. 

Even if we use a computer algebra system, these troubles are serious. (In Ref. [TT], we 
obtain lim n _ oo (0|T fc (M) |0)/(Mn) fc only up to k = 5.) 

To develop an efficient derivation of higher order terms, we pay attention to the fol- 
lowing relations: 

Ft(e) = lim <°Pn°> = Mg> tor , = , (24) 
■— - (Mn) k Q k ' ' ' v ! 



n~ >oo 



where 



/o(6) = sin 2 20, (25) 
g (G) = cos 2 26, (26) 
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f k (Q) = r<W fc -i(6-0) cos 2 20 + gk _ l (Q- ( j ) ) sin 2 20], (27) 

JO 

0fe(6) = /°rf0[^ 1 (e-0)cos 2 20 + / fc _ 1 (0-0)sin 2 20] (28) 
Jo 

for fc = 1,2,.... 

We can prove the above relations from Eqs. (fTKj) and (fTTjl . Both Eqs. (J27)) and (J2%|) 
contain only first-order integrals. Moreover, each of them contains only two terms be- 
ing integrated. Thus, we can compute F o (0), Fi(Q), ... in that order efficiently from 
Eqs. (J2U), (123), (EH), and (J2~%j) using a computer algebra system. (In actual fact, we 
use Mathematica for this derivation.) Eqs. (J27j) and (J25J) constitute a pair of recurrence 
formulas. 

Here, we note some properties of F k (Q). First, /fc(O) and (?fc(0) are analytic at any 
B for k = 0,1,.... In other words, /fc(O) and gk{Q) have Taylor expansions about any 
6o which converge to /fc(O) and <7fc(0) in some neighborhood of Go for k = 0,1,..., 
respectively. We can prove these facts by mathematical induction as follows. To begin 
with, both /o(O) and go(Q) are analytic at any 9 from Eqs. (f23J) and Next, we 

assume that /fc(0) and <7fc(0) are analytic at any 6 for some k € {0,1,...}. Then, 
ifc+i(@) and ^ + i(0) are analytic at any 9 because they are integrals of functions made 
of the sine and cosine functions, /fc(0), and as shown in Eqs. (}2*T|) and (J28j) . Thus, 

by mathematical induction, we conclude that /&(©) and 5^(0) are analytic functions for 
k = 0,1,.... 

From Eq. (J24|) . we can obtain F fc (0) by dividing / fc (6) by fc for fc = 0, 1, .... Thus, it 
is possible that Fk(Q) diverges by marching off to infinity near = 0. However, in fact 
we can show 

FJQ) = = Const.6 2 + 0(9 4 ) for k = 0, 1, (29) 

where Const, denotes some constant. (We prove Eq. (j2Sj) in Appendix 1X1) 

4 Numerical calculations 

In this section, we carry out numerical calculations of (P)(B, x) defined in Eq. (J21)) using 
recurrence relations Eqs. (J27)) and (J28|) . Moreover, we investigate the critical point x c , 
over which the quantum algorithm becomes ineffective for the threshold probability P t h- 

First of all, we need to derive an algebraic representation of (P)(Q,x). We compute 
an explicit form of (P)(Q,x) as follows. First, using recurrence relations Eqs. ()27|) and 
()28|) . we derive fk(Q) and 5^(0). Second, using Eq. we derive Fk(Q) from /^(O). 

Next, using Eq. (J22J), w e derive Cfc(9) from Fjt(O). Finally, using Eq. (|2Tj). we derive 
(P)(Q,x) from Cfc(B), which is the /cth-order term of the perturbative expansion. 

In Ref. jTI], we obtain an explicit form of the matrix element only up to the fifth-order 
perturbation [that is, F 5 (Q)\ because we compute -Ffc(O) from Eq. (fTTj) directly. However, 
in this paper, we succeed in deriving an explicit form of the matrix element up to the 
39th-order perturbation [that is, ^39(0)] with the help of the computer algebra system 
thanks to the recurrence relations Eqs. (|27j) and (f2~SJ). [We do not write down the explicit 
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forms of -F 3 (0), F 39 (9) here except for F 5 (Q), because they are very complicated.] 
By the method explained above, we derive the algebraic form of (P)(Q,x) up to the 
39th-order perturbation. 

However, this explicit form of (P)(9, x) is not suitable for numerical calculation. The 
reason is as follows. Let us consider F^Q) for example. The explicit form of ^(9) is 
given by 

. . 1 45 + 7209 2 - 2569 4 

Fs(0) = 240 + 1 966 0809* ^ 

3 + 3202 + 25604 sm49. (30) 



524 2886 s 

It is very difficult to evaluate the value of F 5 (Q) near B = from Eq. (|30|) directly. If we 
take the limit G — > in the second term of Eq. (J30|) . we obtain 

.. 45 + 7209 2 - 2569 4 

hm — cos49 = +oo. (31) 

e^o 1 966 0809 4 v ; 

However, taking the limit G — > in the third term of Eq. (|30|). we obtain 

.. , 3 + 329 2 + 2569 4 . . ^ 

hm — sm 49 

e-V 524 2889 5 ; 

,. / 3 + 329 2 + 2569\sin49 

= — hm ) 

e^o v 131 0729 4 ; 49 

= -oo. (32) 

As explained above, to evaluate the value of -F 5 (9) near 9 = from Eq. (|371jl directly, 
we have to subtract one huge value from another huge value. Thus, if we carry out 
this operation by computer, an underflow error occurs and we cannot predict a result 
of the numerical calculation at all. [We show that Fk(Q) is analytic at 9 = and 
lim e ^o-^fc(@) — for k = 0, 1, ... in Eq. ([29)1 . However, it is difficult to calculate F 5 (0) 
numerically from Eq. (J30|) .] 

In fact, when 9 = 1.0 x 10" 7 , the second term of Eq. flU is equal to 2.28881 x 10 23 and 
the third term of Eq. (jBUj) is equal to —2.28881 x 10 23 with assuming that the computer 
supports only six significant figures. Hence, the sum of the second term and the third 
term in Eq. (J30|) is equal to zero, and only the first term of Eq. (}30|) . 1/240, contributes 
to ^(9) for 9 = 1.0 x 10~ 7 . However, this numerical calculation is meaningless. We can 
find such an underflow error in almost all the higher order terms, ^(9), ^(9), Fs(Q), 

To avoid this trouble, we carry out the following procedure. We expand the explicit 
form of Cfc(9) in powers of 9 up to the 40th-order term and define (7^(9) as the finite 
power series obtained in the variable 9 for k = 0,1, 39. [Cfc(9) is originally an analytic 
function and it has a Taylor expansion about 9 = 0.] We substitute these {Cfc(9) : k = 
0, 1, 39} for Eq. (J2U) and obtain 

39 i 

(P)(Q,x) = Y,C k (Q)^x k . (33) 
fc=o K - 
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Figure 4: 30th, 40th, and 50th-order polynomials, which we obtain as parts of Taylor 
series of C*4o(0), as functions of for < < (9/5)n. A dashed line, a thick solid line, 
and a thin solid line represent the 30th, 40th, and 50th-order polynomials, respectively. 

We use this (P)(0, x) for numerical calculation. [(P)(0, x) is a polynomial, whose highest 
power in is equal to 40 and whose highest power in x is equal to 39.] 

Here, we make some comments on our approximation method for Cfc(0). In this pa- 
per, we use a polynomial of high degree for approximating (7^(0). The reasons for this 
choice are as follows: (1) to obtain the Taylor series of Cfc(0) is easy, and (2) because 
we can calculate integrals and derivatives of polynomials with ease, Ck{Q) is suitable for 
applying Newton's method. (We use Newton's method for calculating x c later.) However, 
approximation with a polynomial of high degree sometimes causes oscillations, and con- 
sequently errors of numerical calculation. Pade approximant method is effective in the 
treatment of this problem. However, we do not use this method here, because we have to 
carry out tough calculations for deriving the Pade approximants of (7^(0). 

We use a 40th-order polynomial for approximating £7^(0) in this paper. Figure 0] 
shows 30th, 40th, and 50th-order polynomials, which we obtain as parts of Taylor series 
of C*4o(0), as functions of for < < (9/5)tt. A dashed line, a thick solid line, 
and a thin solid line represent the 30th, 40th, and 50th-order polynomials, respectively. 
From Fig. 01 we find that the 30th, 40th, and 50th-order polynomials start to diverge 
near = 3.4, 4.3, and 5.3(radian), respectively. From this observation, we think the 
approximation of C 4O (0) with the 40th-order polynomial, that is C i0 (Q), to be valid in 
the range of < < ir. [Strictly speaking, this is not a rigorous proof but evidence that 
the polynomial expansion up to the 40th-order term is sufficient for approximating Cfc(0) 
for k = 0,1, 39 in the range of < < ir.] 

To investigate the range of x where our perturbative approach is valid, we need to 
estimate the 40th-order perturbation. From numerical calculation, we obtain 



for < < 7i. (From now on, we limit to < < 7r for our analysis, because the 



0< 1-^^40(0)1 < 1.24 x 10 
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Figure 5: (P)(6, x) as a function of x with fixing 9 = max , where max = (l/2)(2M max + 
1)0, M max = 17, # = arcsin ^/l/2 n , and n = 9. Both (P) and x are dimensionless. A thin 
solid curve represents (P)(0,x) obtained by numerical calculation up to the 39th-order 
perturbation. Black circles represent results obtained by Monte Carlo simulations of the 
n = 9 case (nine qubits) with M max = 17. Each circle is obtained for x = 2M max np = 306j9, 
where p is varied from p = 2.0 x 10~ 3 to p = 2.0 x 10~ 2 at intervals of Ap = 2.0 x 10 -3 . 
In these simulations, we make 50 000 trials for taking an average. 

approximation of Cfc(G) with the 40th-order polynomial is reliable in this range, as shown 
in Fig. I3J) Hence, if we limit x to < x < 10.0, the 40th-order perturbation is bounded 
to 



[From now on, we write (P)(Q,x) as the approximate form (P)(Q,x) for convenience as 
far as this naming does not create any confusion.] 

Let us investigate (P)(0, x) obtained in Eq. (}3*3*|) by numerical calculations. To confirm 
reliability of our perturbation theory, we compare the obtained (P)(B,x) with results of 
Monte Carlo simulations of our model in Figs. and El In these simulations, setting 
n = 9 (nine qubits), we fix p and cause phase flip errors (a z errors) at random in each 
trial. We take an average of (0|p (M) |0) p , the probability of observing |0) at the Mth 
step [M = 0, 1, ...,M max (= 17)], with 50 000 trials for each certain value of p. [Because 
(tt/4)V2? = 17.7.., we put M max = 17.] 

Figure shows (P)(@, x) as a function of x with fixing — ©max? where ©max — 
(l/2)(2M max + 1)9, M max = 17, 9 = arcsin ^1/2", and n = 9. (Hence, the independent 
parameter is only p actually.) At x — 0, there is no error in the quantum process and (P) 
is nearly equal to unity. As the error rate x becomes larger, (P) decreases monotonically. 

Figure El shows (P)(@, function of G with fixing p. Because we use the variable 

x = 2Mnp instead of p in the perturbation theory, we have to rewrite x as 



< I— CM)x 

- I 40 , 40V ) 



< 1.24 x 10 
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Figure 6: (P)(G, function of B (radian) with fixing p. Both (P) and B are 



dimensionless. To estimate (P)(0,x), we put x = 20(arcsin y l/2 n ) where n = 9. 
Four thin solid curves represent p = 2.0 x 10 -3 , 4.0 x 10~ 3 , 6.0 x 10~ 3 , and 8.0 x 10 -3 in order 
from top to bottom. Black circles represent results obtained by Monte Carlo simulations 
of the n = 9 cases (nine qubits). Each circle is obtained for = (1/2)(2M + 1)6, where 
9 = (arcsin yl/2™)- 1 , n = 9, and M e {0, 1, M max (= 17)}. 

which we obtain by substituting = lim n _ +00 M9 for x = 2Mnp without taking the limit 
n — > oo, and we give some finite n to Eq. (jHfijl . In Fig. |H1 we set n = 9 and plot curves 
with p = 2.0 x 10~ 3 , 4.0 x 10~ 3 , 6.0 x 10~ 3 , and 8.0 x 10~ 3 in order from top to bottom. 
We also plot results of the simulations. When we plot a result of the simulation for the 
Mth step, we put 



where n = 9 and M e {0, 1, M max (= 17)}. We obtain Eq. ((3ZJ) from Eq. (fTHJ) and 
= lim™ M9. 

From Fig. |3 we notice that the maximum value of (P) is taken at < 7r/4 for 
each p and the shift becomes larger as p increases. This fact means that 0th(Pc,Pth) 
becomes smaller than tt/4, as P t h decreases. [We write Qth{p, Pth) = lim n ^oo M th (p, P t h)0 
and M t h(p, Pth) represents the least number of the operations iterated for amplifying the 
probability of |0) to P t h under the error rate p.] 

Finally, we compute function of P t h- We show the result in Figs. |21andEl We 

obtain x c for < VP t h < 1 as follows. We calculate t h(x,P t h) for given P t h varying x 
from zero, where t h(a;, P t h) = lim^oo M th (x, P t h)0 and M th (x,P t h) represents the least 
number of the operations to amplify the probability of |0) to P t h under given x. (We use 
Newton's method for obtaining a root of for the equation (P)(0, x) = Pth for given x.) 
When x becomes a certain value, we cannot find a root of 0th(^, Pth) and we regard it as 
x c . By repeating this calculation for many points of Pth, we obtain the curve shown in 
Figs. El and H 




= (1/2)(2M + 1)0 = (1/2)(2M+ 1) (arcsin y^I/2™)- 1 , 



(37) 
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Using Eq. (J2TJ), a tangent at P t h = 1 is given by 



x c = c(l - P t h), 



1 8 



(38) 



c = 



Cx(7r/4) 5' 



because Oth(^c-Pth) = vr/4 and x c = for P t h = 1- This means that the algorithm is 
effective for 2Mnp < (8/5) (f — P t h) near P t h = 1, as shown in Fig. El This result is 
similar to a work obtained by Bernstein and Vazirani ^2], and Preskill [HI], as explained 
in Sec. [T] Moreover, we notice x c > —(8/5) log e P t h near P th = from Fig. El 

5 Discussions 

From Fig. 121 we find that the algorithm is effective for x = 2Mnp < (8/5)(l — P t h) 
near P th = 1, and this relation is applied to a wide range of P t h approximately. Thus, 
if we assume that P t h is equal to a certain value (1/2 < P th < 1, for example), we can 
expect that the algorithm works for x = 2Mnp < 0(1) (x is equal to or less than some 
constant.) Hence, if the error rate p is smaller than an inverse of the number of quantum 
gates (2Mn)~ 1 , the algorithm is reliable. If this observation holds good for other quantum 
algorithms, it can serve as a strong foundation to realize quantum computation. 

After we studied decoherence in Grover's algorithm with a perturbation theory in 
Ref. [H], some other groups have tried similar analyses. Shapira et al. investigate perfor- 
mance of Grover's algorithm under unitary noise They assume the noisy Hadamard 
gate and estimate the success probability to detect a marked state up to the first order per- 
turbation. Hasegawa and Yura consider decoherence in the quantum counting algorithm, 
which is a combination of Grover's algorithm and the quantum Fourier transformation, 
under the depolarizing channel fTS]. 
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To put it more precisely, we can obtain the following relations in which Eq. ()39|) is included: 



A Proof of Eq. (ES) 



In this section, we prove Eq. (}2T?|) . which we can rewrite in the form 
f k (Q) = Const.0 fc+2 + O(0 fe+4 ) for k = 0, 1, .... 



(39) 



fk(0) 



oo 



2>f 

3=0 




(40) 



b { h) e k + b < t ) e k+2 + b^Q k+A + ... 



b (k)Qk+2 j 



(41) 



j=0 

for k = 0,1, 
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where f k (Q) and g k (Q) are defined in Eqs. (J23), 423) > and 

We prove Eqs. (jlUj) and (}4*Tj) by mathematical induction. First, when k = 0, we obtain 
the following results from Eqs. (J25j) and 

/o(8) = ^2BHE^e-f 

= 4 e 2 -^e 4 + ^e 6 + ..., (42) 

3 45 

, (6) = cos 2 2e = E^^e 2 T 

1 6 

= l-49 2 + — 9 4 + .... (43) 
3 

Thus, Eqs. pOj) and ({H} are satisfied for fc = 0. 

Next, assuming that Eqs. (J4*U|) and (j4*Tj) are satisfied for some k, we investigate whether 
or not Eqs. (}4T)j) and (}4~Tj) hold for (k + 1). Let us consider Eq. (J4*U|) for (k + 1). From 
Eq. (J2ZJ), we obtain 

/ fc+1 (6) = [° d<p[f k (Q - 0) cos 2 20 + ^(9 - 0) sin 2 20]. (44) 

Here, we expand cos 2 20 and sin 2 20 as follows: 

oo 

cos 2 20 = X^ 2i > ( 45 ) 



j=0 

oo 



sin 2 



20 = ^^0 2 ^' +1 ). (46) 

3=0 

From Eqs. (gOJ), flU}, (03), and pSJ). we can rewrite Eq. (jUJ) in the form 

„Q OO OO 

We) = / EK (fc) c,(© - 0) fc+2(m) 2j 

Jo i=0 j=0 

+6f ) d i (e-0) A + 2 > 2 ^+ 1 )]. (47) 

Applying the following formula 

d0(6 - 0)> 2J = HfZZ: e i+2j+1 (48) 

to Eq. (HU), we find that /fc+i(6) includes only terms of Q k+3 , 9 fc+5 , 6 fc+r , .... Therefore, 
Eq. (jHI holds for (k + 1). Next, let us consider Eq. (JUJ) for (fc + 1). From Eq. lf2Bjl . we 
obtain 

g k+ i(e) = [ d<j)[g k {Q - 0) cos 2 20 + f k (Q - 0) sin 2 20]. (49) 
jo 



Using Eqs. (jUJ), (jUJ), (J3HJ). and pHjl . we can rewrite Eq. (JUJ) in the form 

„Q OO OO 

^i(e) = / «^EE^(e-W 

J ° i=0 j=0 

+af - 0) fe+2 ( i+ V°' +1) ]- (50) 
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Applying Eq. pSjl to Eq. (jHOJ), we find that g k+ i(Q) includes only terms of fc+1 , 6 fc+3 , 
G fe+5 , .... Therefore, Eq. ()41|) holds for (k + 1). Hence, from mathematical induction, we 
conclude that Eqs. (j4Ti|) and (j4Tj) are satisfied for k = 0, 1, .... This implies that Eq. (J3T?|) 
holds. 
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